Kinetic Monte Carlo simulation of shape transition of strained quantum dots 
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The pyramid-to-dome transition in GeajSii-^; on Si(lOO) initiated by step formation on pyramidal 
quantum dots is atomistically simulated using a multistate lattice model in two-dimensions incorpo- 
rating effective surface reconstructions. Under quasi-equilibrium growth conditions associated with 
low deposition rates, the transition occurs at island size Uc following ^/n^ ~ x~^'^^ independent of 
temperature and deposition rate. The shape transition is found to be an activated process. Results 
are explained by a theory based on simple forms of facet energies and elastic energies estimated using 
a shallow island approximation. An asymptotic scaling relation nj^'' ~ for a; — >■ applicable 
to d = 2 or 3 dimensions is derived. The shape transition energy barrier can be dominated by the 
interface energy between steep and shallow facets. 

PACS numbers: 81.15.Aa, 68.65.Hb, 81.16.Dn, 81.16.Rf 



I. INTRODUCTION 

The self assembly of quantum dots in heteroepitaxy 
exhibits very interesting physics and has possible appli- 
cation to device fabrication [J-Q. Growth of Ge or GeSi 
alloy nanostructures on Si(lOO) is the prototype example 
most widely studied. Under typical deposition condi- 
tions, pyramidal islands bounded by shallow (105) facets 
form spontaneously on a wetting layer. They can then 
undergo transitions into multi-faceted dome islands dom- 
inated by much steeper (113) facets and bounded also by 
other facets This shape transition gives rise to a bi- 
modal island size distribution with enhanced dome size 
uniformity [Bj . An atomic pathway based on step bunch- 
ing on the pyramids has been identified [1, 0] ■ 

In this work, we report atomistic dynamic simulations 
of the pyramid-to-dome transition using a fast kinetic 
Monte Carlo (KMC) approach based on a multistate 
solid-on-solid model in two-dimensions (2D). Extensive 
simulations under a wide range of conditions are per- 
formed and a simple analytical description is presented. 
A scaling relation for the transition island size generaliz- 
able to three dimensions (3D) is investigated. 

Kinetic simulation of a strained film is much more 
challenging computationally than the unstrained case 
because of the long-range nature of elastic interac- 
tions. First principles calculations 04Il| and molecular- 
dynamics simulations [l2l . [isj have provided important 
information on the energetics of the relevant surfaces and 
steps, but are in general limited to the studies of static 
properties at small system sizes. Continuum simulations 
in contrast are instrumental for investigating large scale 
and long time behaviors [ij, [l^ , but lacking atomic dis- 
creteness, nucleation events associated with island for- 
mation and shape transition cannot be naturally sim- 
ulated. KMC simulations based on lattice models are 
hence unique in allowing large scale atomistic studies on 
the dynamics of strained heteroepitaxy. 

Using a ball and spring lattice model for elastic solids, 
Orr et. al [l6| conducted early KMC simulations of 
strained layers in 2D. In the simulations, the elastic en- 



ergy of the system has to be computed repeatedly in or- 
der to simulate the atomic hopping events responsible for 
the morphological evolution. Applying more advanced 
algorithms for the solution of the elastic problem and 
the sampling of atomic hopping events, simulations us- 
ing larg e and moderate system sizes in 2D [IMl and 3D 
[2(]| - [24| respectively became possible. The model was ex- 
tended recently to model (105) facets jlH and atomic in- 
termixing with substrate atoms [1^ . Alternatively, KMC 
simulations can also be performed efficiently using more 
approximate forms of elastic interactions j27H29{ . 

This paper is organized as follows. Section II explains 
our multistate model for elastic solids which can account 
for both a shallow and a steep facet. The KMC sim- 
ulation results are presented in Sec. III. In Sec. IV, is- 
land energies and the island transition rate are calculated 
theoretically. In Sec. VI, a scaling relation between the 
transition island size and the Ge concentration is derived. 
We conclude in Sec. VI with some further discussions. 



II. MULTISTATE SOLID ON SOLID MODEL 

Our model is based on a ball and spring square lat- 
tice model of elastic solids for Gea;Sii_a; on Si JJi]. The 
substrate lattice constant is a. = 2.72 A while the film 
material admits a lattice misfit e = 0.04a;. Each node 
on the lattice represents an atom and it is connected to 
its nearest and next nearest neighbors by elastic springs 
with force constants fc^r = 13.85 eV/a^ and k^N = k^/'i 
respectively. This choice gives the correct modulus cn 
of silicon and a shear modulus constant along tangential 
and diagonal directions. 

In this work, (100), (105) and (113) surfaces must be 
effectively simulated. A leveled surface in the model nat- 
urally accounts for a (100) surface. However, lattice mod- 
els generally lead to islands with a single type of side- 
walls at 45° inclination or steeper |16l - [24| . A multi-state 
extension has been introduced recently in Ref. [1^ to 
effectively model the much shallower (105) facets of a 
pyramid in 2D. We now further generalize it to simulate 
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both shallow and steep facets with slopes 



Sl 



1/5 and §2 = 1/2 



(1) 



comparable to those of realistic facets in pyramids and 
domes. Specifically, atoms are normally represented by 
squares. To effectively model surface reconstructions 
leading to specific facets, we allow local deformation of 
the topmost atoms in the film or substrate into trapezoids 
each characterized by a tilt variable at and an extension 
variable Ki . Here, txi gives the slope of the upper surface 
of an atom and equals 



= 0, ±1/5, 



or 



± 1/2 



(2) 



at a locally undeformed region, a shallow facet or a steep 
facet respectively. Allowing atomistically fiat shallow and 
steep facets further requires additional freedoms of ver- 
tical stretching or compression of the topmost atoms by 
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±2/5 
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This characterizes a total of 15 possible local deformation 
states. All lengths are measured in unit of Og throughout 
this paper. Atomic column i with hi atoms is thus trape- 
zoidal in general with the left and right edges of heights 
/i° and h'^ given by 



hi 
h. 
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(4) 
(5) 



A surface step in between column i and i + 1 has a height 

5^=\h1+,-h\\ (6) 

projected along the vertical direction. Figures [lla)-(b) 
show examples of atomic configurations. 

Misfit induced elastic strain is assumed to be com- 
pletely independent of the local deformations associated 
with faceting introduced above. The elastic relaxation 
energy Es of the system is defined as the total energy 
storied in all springs at mechanical equilibrium compared 
with that in the homogeneously strained state. The bond 
energy of the system is defined relative to that of a fiat 
surface by 

{5^)\+ER (7) 

i 

where the facet-type label cti, depending on | |, indi- 
cates if column i is locally undeformed (a^ = 0) or cor- 
responds to a shallow [ai = 1) or steep facet (a^ — 2). 
The facet formation energy per site 4>ai equals = 0, 
^1=5 meV, or (f)2 = 50 meV. These values control the 
relative stability of the facets in our simulations. They 
are chosen empirically so that shallow and steep facets 
start to emerge on islands of appropriate sizes. Also, the 
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FIG. 1: (a) A shallow facet with steps leading to (b) a 
steep facet from a small-scale simulation. In (b), the first 
6 surface atoms from the left have local deformation states 
(cri,Ki) = (|,§),(^,-f),(f,-^),(|,i),(|,-3) and {\,\)- 
Surface atoms in shallow (steep) facets are shaded in red 
(green), while bulk Ge (Si) atoms are shaded in light (dark) 
blue. 



facet interface energy -0(1, i + 1) is non-zero only at the 
boundary between either different facet types or differ- 
ent facet orientations (i.e. ai ^ <^i+i) where it equals 
i>a,ai+i with Vol = i'll = 0.35 eV, 0^12 = V'22 = 0.5 eV, 
and -002 = "001 + "012, assuming Vaa' = ipa'a- The 
ll) term represents surface step energy. It equals "fSi/2 
on a locally undeformed region where 7 = 0.5 eV is 
the nearest neighboring bond energy. If the site i or 
i + 1 belongs to a shallow or a steep facet, it equals 
/3a[l + %— X6xp(l — 6i/sa)] +j{Si — Sa)/2 whcre a is the 
larger of ai and a^+i. Here, f3i = 0.3 eV and (32 — 0.2 eV 
represent the single step energies on shallow and steep 
facets respectively and x = 0.3 dictates the tendency 
of step bunching. To discourage very steep regions, the 
intrinsic step height defined by (5^ =| /li+i — hi \ disre- 
garding local deformation is constrained to (5^ < 1 and 
furthermore each pair of consecutive upward or down- 
ward intrinsic steps with S'i = 1 contributes 0.15 eV to 
the repulsion energy Eji. The constraint also limits the 
step heights 6i to bounded values, although double steps 
in particular, which have heights 2/5 and 1 on shallow 
and steep facets respectively are still possible. 

The KMC approach simulates every hopping event of 
a topmost film atom m according to the rate 



F(m) — Rq exp 



AEb{m) + AEs{m) + E'„ 



(8) 



Atoms then lands random on any other site at most 8 
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columns away. Here, /S.Eb{m) and AEs{m) denote the 
change in the bond energy Ei, and the strain energy Eg of 
the system when the site is occupied versus unoccupied. 
We put E'^ = -7-0.67 eV and Ro ^ 4.1 x 10"s-i. This 
gives the appropriate adatom diffusion coefficient for sil- 
icon (100). Due to the long-range nature of elastic inter- 
actions, the repeated calculations of AEs{m) dominates 
the simulation run time and we handle it using a Green's 
function method together with a super-particle approach 
[23} . Exposed substrate atoms are not allowed to hop. 
Elastic couplings of adatoms with the rest of the system 
are weak and are neglected for better computational ef- 
ficiency. Atomic hoppings are assumed to preserve the 
local deformation states. After every period r, a set of 
deformation states will be updated. We put r — 2 /Tad 
where Tad is the adatom hopping rate on an locally un- 
deformed region easily calculable from Eq. ([8]). At an 
odd (even) numbered updating event, all odd (even) lat- 
tice sites will be updated. The variables ct^ and Ki at 
those sites are re-sampled from the allowed set of 15 pos- 
sible combinations using a heat bath algorithm based on 
the relative probability exp{—Eb/kT). Our model obeys 
detailed balance. The dynamical rules described above 
reduces back to those used in Refs. [l^l at locally unde- 
formed regions. We have critically checked our software 
implementation, in particular using a Boltzmann's dis- 
tribution test [23], which is found to be indispensable in 
verifying that practically all, but not only the dominating 
hopping pathways can be correctly simulated. We also 
have checked in small scale simulations that restricting 
hoppings to only nearest neighboring sites rather than 
allowing long jumps gives similar results except for an 
insignificant shift in the time scale. Wetting layers on 
the substrates are believed to be relatively immobile and 
are not simulated for simplicity. 



III. SIMULATION RESULTS 

Figures [Ija) show a shallow facet with steps from 
a small scale simulation. The steps are subsequently 
smoothed out by the formation of a steep facet as shown 
in Fig. [Ijb) . Figure [2][a) shows snapshots of a surface 
from a typical simulation of deposition at 4 ML/s, 700°C 
and a; = 1 on a substrate of width 2048. Successive pro- 
files are displaced vertically for clarity. Stepped mounds 
first develop and some of them matures into pyramids 
bounded by shallow facets as explained in Ref. [1^. 
Some of the pyramids further turns into domes bounded 
mainly by steep facets often with regions of shallow facets 
at the top. Figure [2](b) shows the detailed evolution of 
one of the domes. Steep facets on either side of a pyra- 
mid form independently. The transition hence often goes 
through a meta-stable half-dome stage. The formation of 
most domes is preceded by steps appearing on the shallow 
facets as shown in Fig. [T]and Fig. [Hb) and as proposed 
in Ref. Q. A close examination reveal that these steps 
are highly dynamic and continuously bunches, separate 
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FIG. 2: (a) Surface profiles at 1 to 6 ML coverage simu- 
lated at 700° C with deposition rate 4 ML/s on a lattice of 
width 2048. (b) Detailed profiles at 4 to 4.5 ML coverage 
showing a pyramid-to-dome transition corresponding to the 
leftmost dome in (a). In (a) and (b), each successive surface 
corresponds to the deposition of a further 1/4 and 1/64 lay- 
ers respectively and is displaced vertically for clarity. Locally 
deformed regions, shallow facets and steep facets are shaded 
in blue, red and green respectively. 

and diffuses around. After accumulating a considerable 
total step height, they transform highly reversibly into 
a steep facet. As the total height of the steps increases, 
the resulting steep facet is more stable and eventually 
become fully stabilized. A smaller number of domes are 
initiated instead by the formation of steep facets at the 
base of the pyramid when shallow facets temporarily de- 
cay into unfaceted regions due to thermal excitations. 

Large scale deposition simulations have been per- 
formed at temperature T from 450° C to 850° C at a; = 1 
on lattices of width 2048. The deposition rate R varies 
from 0.006 to 8 ML/s chosen empirically to generate typ- 
ically 3 to 5 pyramids or domes on each substrate. The 
low island density minimizes elastic interactions among 
islands which are known to alter the dynamics [13, |3l| . 
An island is defined as one in which all constituent 
columns must be at least 2 atoms tall. We measure island 
size in number of atoms n so that ^/n is proportional to 
the linear size of the island. Also, island aspect ratio is 
defined by r = h/2l where h is the height of the high- 
est point of the island and 21 is its basewidth. Figure 
[3la) shows a scatter plot of the aspect ratio r against 
^/n for all islands from 3 independent runs at each tem- 
perature. Measurements are conducted throughout the 
evolution. Time averaging of values associated with in- 
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FIG. 3: (a) Plot of island aspect ratio r vs root of island size 
■^n for various T and R from simulations of deposition with 
lattice width 2048. (b) Plot of r vs -Jnjx'^ with C = 1.69 
for various x from similar simulations. Inset: Log-log plot of 
^Jri^ vs X where ric is the transition island size. The dotted 
line shows a linear fit to the data giving ^ = 1.69. The solid 
line represents a theoretical result. 



dividual islands over short periods are performed, but 
no ensemble averaging is done as each island develops 
in general at a different pace. We observe that r first 
converges towards 0.1 as islands transform from stepped 
mounds into pyramids. It then rises again to around 0.2 
characterizing the dome transition similar to experiment 
findings in Ref. . The much lower density of the data 
points at 0.1 < r < 0.2 corresponding to highly unsta- 
ble intermediate states was also observed in Ref. 0. 
The morphologies of these intermediate states have been 
shown in Fig. 2(b). Results in Fig. 3 reveals two distinct 
trends. For the runs at T > 650°C, all islands follows an 
identical evolution path with the dome transition occur- 
ring at size Uc — 900. In contrast, at lower temperature 
T < 600° C, the transitions are delayed randomly to in- 
creasingly larger sizes. At T = 450° for instance, ric 
ranges from about 900 to 2000. We will explain these 
distinct trends at the end of this section. 

Similar deposition simulations are also performed at 
T = 700° C and i? = 1 ML/s for Ge concentration x 
from 0.6 to 1 with 3 independent runs in each case. We 
find that the dome transitions occur at increasingly larger 
island sizes as x decreases. The precise dependence is 
easily illustrated by a rescaled plot of r against ^Jnlx~'-' 



with C = 1.69 as shown in Fig. [3jb). Data for various 
values of x collapse reasonably well into a single curve, 
implying 



J{n/x 



(9) 



where / is a rescaled function. To calculate C, used above, 
we have measured the transition size ric by averaging the 
sizes of all islands right at the transitions with 0.12 < 
r < 0.16. The resulting plot of against x in log-log 
scales is shown in the inset of Fig. Elb). A linear relation 
observed in the log-log plot implies 



(10) 



and a linear fit gives C = 
be explained in Sec. IV. 



1.69. This scaling relation will 



pyramid 

half-dome ^ 

dome X ^y-{:' 




0.0015 



0.1 0.12 0.14 0.16 0.18 0.2 0.22 
Aspect Ratio 

FIG. 4: Plot of bond energy Eb, elastic energy Eg and to- 
tal energy E against island aspect ratio r from the annealing 
of an initially pyramidal island on a substrate of width 256 
at 550° C (symbols). The solid lines show theoretical results. 
The schematic diagram shows a pyramid with additional lay- 
ers on one of the shallow facets during the transition into a 
half-dome. Inset: An Arrhenius plot of the dome transition 
time r. 

To study the energies of individual islands, we have 
performed simulations on annealing of single pyramids 
directly constructed on smaller substrates each of width 
256. The annealing temperature is 550°C. The pyramid 
is initially bounded by shallow facets and sits on an oth- 
erwise empty substrate surface. It contains 1230 atoms 
and has a basewidth slightly less than the lattice width. 
This number is chosen empirically to be just sufficient to 
ensure an irreversible dome transition. Figure |4] shows a 
scatter plot of the system bond energy E'f,, strain energy 
Es and total energy E — Ef, + against r measured 
during the annealing from 16 independent runs. Since a 
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single island dominates, these energy of the whole system 
approximates those of an island. Only time averaging 
of the values over short periods but no ensemble aver- 
aging has been carried out. The symbols used indicate 
if the data points correspond to pyramids, half-domes, 
or domes. The geometries are identified reliably by the 
number of steep facets present. The result indicates that 
there is an energy barrier for the transition. Moreover, a 
number of data points associated with half-domes clus- 
ter around r ~ 0.12 — 0.14 showing that the geometry 
characterizes a meta-stable state. 

We next show that the dome transition is an activated 
process. We have repeated the above simulations on the 
annealing of pyramids at T from 450°C to 850°C. A dome 
transition time r defined as the average annealing dura- 
tion required to reach an aspect ratio r > 0.12 is mea- 
sured. Values of t each averaged over 16 independent 
runs are plotted against l/T in the inset in Fig. 01 The 
data fits well to 

T = Toexp{no/kT) (11) 

with rJo = 1-97 eV and To = 2.9 x 10"^^ s. The Arrhenius 
temperature dependence is typical of activated processes. 
The value of Oq will be explained in Sec. IV. We further 
repeat the simulation 300 times at 700° C. The values of r 
measured are histogrammed. The result is well fitted by 
an exponential distribution. This further supports that 
the dome transition is an activated process. 

With the shape transition time r known, the distinct 
trends followed by the data in Fig. ^a) can now be ex- 
plained. The deposition rate R has been empirically cho- 
sen to produce a constant and sparse island density. The 
choice hence ultimately depends on pyramid nucleation 
and coarsening dynamics and is not directly related to 
the dome transition dynamics. For the T > 650°C runs, 
we find that 1/r ^ R. The shape transition is thus 
fast compared with deposition and hence also the island 
growth. The dome transition process is only limited by 
the availability of atoms. The geometry as characterized 
by the aspect ratio r therefore only depends on i/n and 
is independent of T and R as shown in Fig. ^a) . In con- 
trast, for the T < 650°C runs, we get l/r < R. Island 
growth can then be fast enough to out-run the dome tran- 
sition, which becomes rate-limited. There is a significant 
random waiting time for the dome transition following 
an exponential distribution during which the island may 
already have grown to a larger size. The transition thus 
occurs at a more broadly distributed island size ric as 
observable from Fig. ^a). Note that if we consider for 
example a higher island density by increasing the values 
of R used, the characteristic temperature separating the 
two trends, which is found to be around 650°C here, will 
increase. 

IV. THEORY OF SHAPE TRANSITION 

We now present a detailed theoretical analysis based 
on generic forms of elastic and facet energy for the transi- 



tion of a pyramid into a half-dome in 2D, which is applied 
to interpret our KMC simulation results. Half-domes are 
meta-stable and they quickly transform into domes. Our 
formulation is consistent with that of Montalenti et al. 
who have shown using energy parameters for Ge/Si from 
first-principle calculations that the dome transition is en- 
ergetically favorable for sufficiently large pyramids [6]. 

Consider an island of size n initially in the form of a 
pyramid with a half-basewidth Iq. We have 

n = sill. (12) 

Geometrical rearrangements can lead to the formation 
new atomic layers of total vertical thickness v on one of 
the facets as shown in the schematic diagram in Fig. 2] 
The new half-basewidth I is then related to v by 

n = siP + [ub — ua)v (13) 

where ua and u b denote the positions of the midpoints 
A and B on the edges of the new layers measured from 
the apex of the base pyramid. A single atomic step on a 
shallow facet has a height si. We assume for simplicity 
that there are v/si single steps at point B. Using Eq. ([7]), 
the bond energy of the pyramid is 

£;f = 201/ -f2i/;oi -1-^11 (14) 

The pyramid becomes a half-dome when all the steps at 
point B turn into a steep facet of A = v/{s2 — si) columns 
wide. The bond energy of the resulting half-dome also 
follows from Eq. ^ and we get 

j^hdome ^ 2cj),l + 2^Poi + i'll + 't^^v + 2V'i2 (15) 

S2 - Si 

Neglecting the small difference in the elastic energies of 
the two geometries, the island takes the form with the 
lowest bond energy Ef, — minimum{E^^ , E^'^°™^} . 

The elastic energy of the pyramid and the half-dome is 
assumed to be identical and is calculated by approximat- 
ing both edges of the new layers as vertical walls located 
at A and B. A shallow island approximation (32l.[33j gives 

Es=Ce^J J dxdx' s{x)s{x') In * ~ * , (16) 

where s{x) denotes the local surface slope of the is- 
land at position x and Cc is a spatial cutoff. We put 
s{x) — ~sgn{x)si+vS{x~UA)~'v6{x'~UB) where sgn{x) 
and S{x) represent the sign function and the Dirac delta 
function. In 2D, C = ala^/ire^Y where ct;, cx e is the xx 
component of the bulk misfit stress and Y is the Young's 
modulus. For our model, a simple calculation based on 
lattice elasticity gives C — 4A;„ag/37r. 

Performing the integrations in Eq. (llOp . we have 
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2;2 



p=A,B 



Hn 



y2 In 



ub - 



(17) 



where = — 1 and = 1. In the fohowmg, we put 
Oc = e~'^/^A where A = (A + 2ua)/'2. is the average spa- 
tial extent of the misfit force monopoles apphed over 
the edges at A and B. It can be shown that this choice 
gives the correct energy when approaching the point force 
hmit. 

From simple geometry, ua — v/Asi. We calculate ub 
by minimizing the total island energy E at small v. Lin- 
earizing E from Eqs. P^ - (IT5|) and (|17p w.r.t. v, it can 
be shown after some algebra that for both pyramids and 
half-domes, E is minimized at 



Ub = h 



1 + 4 exp - 



01 



-1/2 



(18) 



The total energy cost A£' of an island compared with the 
initial pyramid can then be calculated. For a half-dome, 
we get, up to linear order in v, 



AE = 



32 - Si 



2Ce^silo In 



Iq + Ub 

lo — UB 



2tP 



12- 



(19) 

Equation (fT8|) gives the energetically most favorable po- 
sition for the initial formation of a steep facet. Using the 
KMC model parameters for a; = 1, it gives ub/Io — 0.53 
for islands around the transition size. The result is in 
general close to a limiting value ub/Iq = 5~^/^ ~ 0.447 
obtained by neglecting the shallow facet formation en- 
ergy 01 . 

The energies Et, Eg and E are numerically calculated 
for various layer thickness v adopting the KMC model 
parameters for a; = 1 using Eqs. (HSll-lIISl) and (fT7l) - (fT8l) . 
The island aspect ratio r is also calculated as a function 
of V using r = si/2 + v/Al and Eq. In Fig. IH 

the energies are plotted as solid lines against r. We have 
assumed an island size n = 1183 atoms measured during 
the dome transition in the KMC annealing simulations 
responsible for the data points in Fig. The only tun- 
able parameter is a fitted additive constant 3.4 eV for 
Eb, which accounts for the bond energies of all other ex- 
citations in the system. It nevertheless plays no role in 
further calculations. The theoretical estimates of the en- 
ergies generally show reasonable quantitative agreement 
with simulation results as observed in Fig. 21 The main 
discrepancies are due to errors in Es , since the shallow is- 
land approximation is known to overestimate the elastic 
relaxation at large r. Nevertheless, important features 
including a shape transition energy barrier and a meta- 
stable half-dome state are correctly reproduced and these 
will be further studied. 

From Fig. SI both theory and KMC simulation show 
an energy barrier for the dome transition associated with 



a maximum in the total energy E. Its location follows 
theoretically from i?™ = E^'^°"^'^. We get a barrier height 



AE„ 



0.88 eV which occurs at r = 1.03 or v = 0.81. 



It corresponds to v/si ~ 4 new atomic layers on the 
shallow facet. For r below and above 1.03 respectively, 
pyramid and half-dome are the energetically preferred 
states. Due to the small value of v at the barrier, AEmax 
is dominated by the steep facet interface energy term as 
can be deduced from Eq. ([T^ . i.e. AEmax — 2V'i2 = 
1 eV. The dominance of the steep facet interface energy 
on the transition energy barrier may be a general feature 
applicable also in 3D. 

The existence of an energy barrier confirms that the 
dome transition is an activated process as suggested in 
Sec. III. The transition rate R hence follows the Arrhe- 
nius form R — v e:xjp{~' AEmax / kT) , where v denotes an 
attempt frequency. Assuming that the transition is lim- 
ited by the diffusion of adatoms on the shallow facet, 
one expect v cx pD, where p = exp{—Ead/kT) and 
D cx exp(— 51arf//cT) are the adatom density and diffu- 
sion coefficient on the shallow facet. Here, Ead = 0.6 eV 
and ~ 0.57 eV are the adatom formation energy and 
hopping energy barrier on the shallow facet calculated 
using Eqs. ^ and ([8]). In particular, fiad is not far from 
a previous estimate from first-principles calculations (sj] . 
The dome transition time r cx 1/i? is hence given by 



T cx exp 



AEr, 



Ead + ^ac 



kT 



(20) 



A comparison with Eq. (ITU) leads to fto — AEmax + 
Ead + ^ad- It givcs flo — 2.05 cV in agreement with 
51o = 1-97 eV obtained above from KMC simulations. 



V. SCALING OF SHAPE TRANSITION SIZE 

We first assume quasi-equilibrium conditions in which 
the dome transition is fast compared with island growth. 
It applies to the cases of slow deposition and small tran- 
sition barrier and is valid for our KMC simulations at 
T > 650° C (see Fig. E^a)). The island energy E from 
theory as shown in Fig. 2] exhibits a local minimum rep- 
resenting the meta-stable half-dome state. The energy 
rises again at larger r because the base pyramid then 
becomes too small to relief the elastic energy efficiently. 
The dome transition occurs only if the half-dome is suffi- 
ciently stable, say of energy kT below that of the initial 
pyramid. For island at the transition size ric, the mini- 
mum of AE hence follows AEmin = —kT. We can then 
calculate Uc numerically using Eqs. (|13p . ((T5|) . p7)) and 
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(fT8| and the values are plotted against solid line 

in the inset of Fig. ^h) . Note that no tunable parame- 
ter is involved. The values are in reasonable agreement 
with the KMC results and supports the scaling relation 
in Eq. ([T0|) with C, — 1.49 consistent with 1.69 found in 
simulations. 

In addition to the numerical estimate of the exponent 
above, better insights are obtained by deriving an exact 
exponent C = 2 valid asymptotically in the small misfit 
limit, i.e. e oc x — >■ 0. Assume that the relative position 
of the steep facet is independent of e so that ub oc lo, 
which will be justified later. Simple scaling properties in 
2D elasticity imply that the change in the elastic energy 
of a half-dome compared with the initial pyramid follows 
AEs — —e^lQg2{v/lo) for some function 92- This is also 
explicitly derivable from Eq. p7)) . The total energy cost 
of a half-dome is hence 

AE^A2V~e^llg2{v/lo) + B2 (21) 

where A2 = {4>2 — 4>i) / {,82 — si) and B2 — 21^12. At 
island transition size Uc and considering a layer height v 
minimizing AE to AEmin = —kT, we have 

eHlg2{v/lo)-A2V = B2 + kT (22) 

It means that the elastic energy gain must overcompen- 
sate the facet formation energy cost by an excess amount 
B2 + kT. As e — >■ 0, we will see in the following that 
the shape transition occurs at a larger island size. Both 
energy terms on the l.h.s. of Eq. ([22l) increase unbound- 
edly and must balance each other, while the constant 
energy excess becomes negligible. Therefore, we have 
e^^o52(w//o) — Av. The meta-stable half-dome state at 
transition is thus characterized by the scaling solution 
lo ^ V ^ e^^. Using y/n^ cx Zo and e = 0.04a;, we get 

y/TTc'^x-'^ (23) 

i.e. C = 2. This solution is consistent with the assump- 
tion that Ub is independent of e as deduced using Eq. 
(fT8| . It also justifies that both terms on the l.h.s. of Eq. 
increase unboundedly as e — > 0. 

For finite e and x, the energy excess B2 -t- kT in Eq. 
((22)) is not negligible. It gives a finite size correction to 
the exact scaling in Eq. (f23| . This results at an effective 
scaling y/rIZ ^ x~'^ with C ^ 2 consistent with C, = 1.69 
found in our KMC simulation. 

In 3D, the dome transition is initiated by step bunching 
close to the mid-level on a pyramid [6] . Generalizing our 
discussion above, we consider the formation of a square 
ring of steep facet of vertical thickness u on a 3D pyramid 
of basewidth 2Zo. Generalizing Eq. (PT|). the energy cost 
is 

AE = A:ivlo - e^llg3{v/lo) + B3I0 (24) 

where the terms on the r.h.s. similarly denote the facet 
formation energy, the elastic energy gain, and the steep 



facet interface energy respectively, for some smooth func- 
tion f/3 and constants A^ and B^ independent of v and 

lo- A similar calculation leads to nl^^ oc ~ x^^. This 
generalizes Eq. to 

n,!/'^ - x-( (25) 

with the same exponent C, = 2 in dimension d = 2 or 3 
for x ^ 0. For finite x, an effective exponent C 2 is 
expected in both dimensions. 

The asymptotic scalings derived above for quasi- 
equilibrium conditions essentially follows from the bal- 
ance between the steep facet formation energy and the 
elastic energy gain which scales with the island size dif- 
ferently. It is analogous to the scaling predicted for island 
formation size based on the Asaro-Tiller-Grinfeld insta- 
bility theory [l|. Due to the simplicity, it is also rather 
robust as will be shown below. 

Instead of assuming quasi-equilibrium conditions, the 
dome transition can be limited by the kinetics. This may 
be appropriate in particular at 3D for small x since the 
barrier predicted above can become too large to over- 
come. The transition is then delayed to a larger island 
size which lowers the barrier. Let us then assume a very 
simple transition criterion that the energy barrier AE^ax 
must not exceed a given value, say a few times of kT. A 
similar calculation for a; — )■ again arrives at Eq. (j25p 
with the same exponent C = 2 in both 2D and 3D. The 
solution also requires lo ~ but w ~ and we have 
used gd{z) cx z for z — > which readily follows in the 2D 
case from Eq. ([T^ after neglecting a logarithmic factor. 

VI. DISCUSSIONS 

In the calculations above, x denotes the actual Ge con- 
centration in the film so that e = 0.04j:. In experiments, 
it can differ greatly from the nominal concentration due 
to intermixing with substrate atoms and this complicates 
interpretation of experiment results. As a very rough 
estimate, experiments on deposition of pure Ge at e.g. 
450°C and 700°C have found dome transition occurring 
at island volume 2800 nm^ ^] and 2 x 10^ nm^ [36| . 
Neglecting compositional non-uniformity, it has been es- 
timated that the actual Ge concentrate in the islands is 
X = 0.82 and 0.43 respectively at 450°C and 700°C [13 ■ 
This gives a very preliminary estimate for the scaling 
exponent C ~ 2.2 consistent with the asymptotic value 
C = 2 derived above, although more experiments are re- 
quired for a reliable conclusion. 

Our simulations and theoretical calculations show 
that occurrence of well-defined dome transitions depends 
strongly in particular on the formation and interface 
energies of steep facets. The detailed dependences of 
these energies on the Ge concentration and temperature 
are neglected. We have also neglected the spatial non- 
uniformity of Ge concentration, surface stress and real- 
istic crystal elastic anisotropy. They should have signif- 
icant quantitative impacts on the shape transition, but 
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arc not expected to alter the dynamics described here 
qualitatively. 

In summary, we have generalized a multi-state lattice 
model for elastic solids to account for both shallow and 
steep facets with tunable energy parameters. Using this 
model, we perform kinetic Monte Carlo simulations to 
study the pyramid-to-dome transition in the heteroepi- 
taxy growth of GexSii-a, on Si in 2D. For sufficiently 
slow deposition, the shape transition occurs at an island 
size independent of temperature and deposition rate. A 
scaling relation between the transition size and the Ge 
concentration is observed. For fast deposition, the tran- 
sition can be delayed randomly to a larger island size. For 
annealing simulations, the shape transition time is found 
to follow an Arrhenius form. A theory based on elastic 



energy in the shallow island approximation and simple 
forms of facet formation energies is derived. Numerical 
solutions of the energetic equations give island energies, 
shape transition size, and shape transition rate in reason- 
able agreement with simulations. The shape transition 
energy barrier is dominated by the interface energy be- 
tween the shallow and the steep facets. We have also 
derived analytically an exact scaling rule between the 
transition size and the Gc concentration applicable in 
the limit of small Ge concentration which is expected to 
be valid in both 2D and 3D. A finite size correction to 
the scaling at higher Ge concentration is explained. 

This work was supported by HK GRF, Grant No. 
PolyU-5009/06P and PolyU Grant No. G-U354. 
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